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Abstract 

A  comprehensive  scheme  is  described  to  construct  rational  solid  T- splines  from  boundary  triangulations  with  arbitrary  topology.  To  extract 
the  topology  of  the  input  geometry,  we  first  compute  a  smooth  harmonic  scalar  field  defined  over  the  mesh  and  saddle  points  are  extracted  to 
determine  the  topology.  By  dealing  with  the  saddle  points,  a  polycube  whose  topology  is  equivalent  to  the  input  geometry  is  built  and  it  serves 
as  the  parametric  domain  for  the  solid  T-spline.  A  polycube  mapping  is  then  used  to  build  a  one-to-one  correspondence  between  the  input 
triangulation  and  the  polycube  boundary.  After  that,  we  choose  the  deformed  octree  subdivision  of  the  polycube  as  the  initial  T-mesh,  and 
make  it  valid  through  pillowing,  quality  improvement  and  applying  templates  to  handle  extraordinary  nodes  and  partial  extraordinary  nodes. 
The  obtained  T-spline  is  C2 -continuous  everywhere  over  the  boundary  surface  except  for  the  local  region  surrounding  poly  cube  corner  nodes. 
The  efficiency  and  robustness  of  the  presented  technique  are  demonstrated  with  several  applications  in  isogeometric  analysis. 
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1.  Introduction 

For  tight  integration  of  Design-Through- Analysis,  isogeo¬ 
metric  analysis  [11]  was  proposed  which  utilizes  spline  func¬ 
tions  as  a  basis.  The  current  state-of-the-art  in  engineering  de¬ 
sign  and  analysis  is  built  on  disparate  geometric  foundations. 
Spline  representation  is  popular  in  design  while  polygonal  mesh 
representation  is  generally  used  in  analysis.  This  leads  to  many 
translational  difficulties  which  affect  the  efficiency  and  accuracy 
of  the  entire  process.  Isogeometric  analysis  utilizes  the  same  ba¬ 
sis  functions  as  the  geometry  representation,  and  consequently, 
analysis  can  be  carried  out  over  an  exact  spline  representa¬ 
tion  of  the  geometry.  Similar  to  other  physically-based  anal¬ 
ysis,  solid  models,  which  can  represent  both  boundary  shape 
and  interior  volume,  are  required  for  many  applications  in  iso¬ 
geometric  analysis.  A  fundamental  step  for  the  unified  Design- 
Through-Analysis  technologies  is  to  automatically  construct 
solid  trivariate  spline  models  from  boundary  representations. 

The  advent  of  T-splines  [16]  gives  more  flexibility  for  ge¬ 
ometric  modeling,  allowing  local  refinement,  non-rectangular 
domains  in  2D  and  non-cubic  domains  in  3D.  T-splines  can 
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represent  a  complicated  design  with  complex  topology  as  a  sin¬ 
gle  watertight  geometry,  avoiding  splitting  the  model  into  sev¬ 
eral  patches.  T-splines  are  a  superior  alternative  to  and  also  are 
compatible  with  NURBS,  which  is  the  current  geometry  stan¬ 
dard  in  CAD  systems.  The  flexibility  of  T-splines  makes  them 
an  ideal  discretization  technology  for  isogeometric  analysis. 

Several  methods  have  been  developed  recently  to  con¬ 
struct  solid  T-splines.  A  solid  T-spline  generation  method  was 
described  for  genus-zero  geometries  [4].  In  [13],  trivariate 
T-splines  were  defined  based  on  the  generalized  poly  cube 
parametrization.  Another  spline  scheme  based  on  poly  cubes, 
called  restricted  trivariate  poly  cube  splines,  was  developed 
in  [21].  This  algorithm  is  based  on  semi- standard  T-splines. 
It  requires  calculation  of  weights  and  the  obtained  T-spline 
elements  are  uniform.  For  all  these  methods,  the  polycube  is 
generated  manually  for  a  given  geometry.  In  addition,  the  con¬ 
structed  T-spline  model  may  contain  some  negative  Jacobian 
elements,  which  is  unsuitable  for  analysis. 

In  our  earlier  work,  we  developed  a  mapping-based  ratio¬ 
nal  solid  T-spline  construction  method  for  genus-zero  geom¬ 
etry  from  the  boundary  surface  triangulation  [23].  To  extend 
this  algorithm  to  more  general  geometry,  we  first  extract  the 
topology  of  the  input  geometry  and  build  a  polycube  which  ap¬ 
proximates  the  input  geometry,  by  computing  a  harmonic  field 
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and  dealing  with  its  critical  points.  Due  to  its  regular  structure, 
the  polycube  is  suitable  for  serving  as  the  parametric  domain 
of  the  tensor-product  spline  representations.  Here,  a  trivariate 
solid  T- spline  is  built  upon  the  generated  poly  cube.  We  then 
build  a  parametric  mapping  between  the  triangulation  and  the 
boundary  of  the  poly  cube.  In  the  following  steps,  an  octree  sub¬ 
division  is  applied  to  the  polycube  and  the  initial  T-mesh  is  a 
deformation  of  its  subdivision.  The  subdivision  continues  until 
the  surface  approximation  error  is  less  than  a  given  tolerance. 
After  that,  the  valid  T-mesh  is  obtained  through  pillowing,  T- 
mesh  quality  improvement,  and  applying  templates  to  handle 
extraordinary  nodes  and  partial  extraordinary  nodes.  Finally, 
Bezier  elements  are  extracted  with  all  positive  Jacobians,  and 
they  are  suitable  for  isogeometric  analysis. 

Compared  with  other  existing  methods,  our  solid  T-spline 
construction  scheme  has  several  attractive  properties:  (1)  it  can 
be  used  for  arbitrary  topology  and  the  polycube  is  created  au¬ 
tomatically;  (2)  the  obtained  trivariate  T-spline  has  a  one-piece 
representation  and  it  contains  very  few  irregular  nodes  where 
the  continuity  degenerates;  (3)  it  employs  rational  T-spline  ba¬ 
sis,  which  guarantee  partition  of  unity  by  definition;  (4)  it  pro¬ 
duces  high  quality  analysis- suitable  T-spline  elements  with  all 
positive  Jacobians;  and  (5)  with  an  adaptive  refinement  scheme, 
the  resulting  T-spline  model  is  an  efficient  representation  for 
analysis. 

The  remainder  of  this  paper  is  organized  as  follows.  Section 
2  reviews  related  work.  Section  3  presents  an  overview  of  the 
construction  algorithm.  The  polycube  construction  algorithm 
is  presented  in  Section  4.  The  T-mesh  construction  algorithm 
from  the  polycube  is  described  in  Section  5  and  solid  T-spline 
construction  is  presented  in  Section  6.  Section  7  presents  ex¬ 
amples  and  Section  8  draws  conclusions. 

2.  Previous  Work 

Harmonic  Fields.  A  harmonic  field  is  the  solution  to  the 
Laplace  equation  with  given  boundary  conditions  and  it  defines 
a  scalar  or  vector  field  over  the  domain.  For  a  given  mesh,  har¬ 
monic  fields  can  be  computed  by  solving  a  linear  system  of  alge¬ 
braic  equations  with  imposed  boundary  conditions.  Harmonic 
fields  have  certain  desirable  properties,  such  as  smoothness, 
and  are  free  of  extraneous  critical  points.  Due  to  these  proper¬ 
ties,  harmonic  fields  have  been  shown  to  provide  effective  tool 
for  a  number  of  geometry  processing  problems.  Dong  et  al.  [3] 
traced  the  integral  lines  through  the  gradient  and  orthogonal 
vector  fields  of  a  harmonic  field  for  quadrilateral  remeshing  of 
arbitrary  manifolds.  Based  on  a  harmonic  map,  a  3D  geomet¬ 
ric  metamorphosis  method  was  developed  for  any  two  objects 
which  are  topologically  equivalent  to  a  sphere  or  a  disk  [12]. 
Joshi  et  al.  [19]  utilized  harmonic  coordinates,  which  are  gen¬ 
eralized  bary centric  coordinates,  in  volume  deformation. 

Surface  Parameterization.  A  surface  parameterization  is  a 
one-to-one  mapping  from  one  surface  in  3D  to  a  suitable  planar 
domain.  Parameterization  is  a  powerful  tool  and  necessary  for 
many  geometry  processing  tasks,  including  data  fitting,  texture 
mapping,  and  remeshing.  Many  significant  advances  have  been 
made  for  surface  parameterization  [5,  6, 18].  In  [8],  a  conformal 


mapping  method  was  presented  to  map  a  genus-zero  closed 
surface  onto  the  unit  sphere  by  minimizing  the  harmonic  energy 
of  the  map.  For  parameterization  of  an  arbitrary  genus  object 
to  simpler  surfaces  of  the  same  genus,  the  mesh  is  usually 
first  segmented  into  disk-like  patches  and  then  each  patch  is 
mapped  onto  the  corresponding  plane.  In  [7],  Gu  et  al.  solved 
the  problem  of  global  conformal  parameterization  for  surfaces 
of  arbitrary  topology,  with  or  without  boundaries. 

Polycube  Generation  and  Application.  A  polyube  is  a  solid 
composed  of  cubes.  It  can  be  used  to  very  roughly  approxi¬ 
mate  the  geometry  of  a  3D  object  while  faithfully  replicating  its 
topology.  Due  to  its  highly  regular  structure,  the  polycube  can 
be  used  as  the  parametric  domain  for  surface  parameterization 
and  spline  modeling.  However,  in  practice  due  to  the  complexity 
of  shapes,  poly  cubes  are  usually  constructed  manually,  entail¬ 
ing  considerable  effort.  In  order  to  produce  polycubes  with  less 
user  intervention,  He  et  al.  [9]  developed  a  method  to  construct 
a  3D  polycube  by  extruding  the  axis-aligned  polygons  which 
approximate  the  horizontal  curved  intersection  contours.  Based 
on  the  polycube,  a  global  parameterization  technique,  polycube 
map  [19],  was  first  used  for  seamless  texture  mapping.  Wang 
et  al.  developed  a  technique  to  build  the  polycube  splines  upon 
the  polycube  map  for  surface  modeling  [20].  In  [13],  an  algo¬ 
rithm  was  developed  to  construct  trivariate  T-splines  over  gen¬ 
eralized  poly  cubes  with  a  global  “one-piece”  representation  for 
general  volumetric  data.  In  [21],  a  theoretical  volumetric  mod¬ 
eling  framework  was  presented  to  construct  restricted  trivariate 
poly  cube  splines,  in  which  the  blending  functions  are  strictly 
bounded  within  the  solid  polycube  domain. 

Solid  Spline  Modeling.  Only  a  few  works  have  been  de¬ 
voted  to  solid  spline  modeling.  In  [10],  a  trivariate  simplex 
spline  modeling  method  was  developed  based  on  a  tetrahedral 
decomposition  of  any  3D  domain  with  complicated  geometry 
and  arbitrary  topology.  In  [24],  a  skeleton-based  method  was 
developed  to  construct  solid  NURBS  for  isogeometric  analysis 
of  arterial  blood  flow.  A  method  was  presented  in  [1]  to  gener¬ 
ate  NURBS  parameterizations  of  swept  volumes  by  sweeping  a 
closed  curve,  and  isogeometric  analysis  was  applied  to  the  gen¬ 
erated  NURBS  model.  Based  on  discrete  harmonic  functions, 
a  volumetric  parameterization  was  used  to  construct  a  single 
trivariate  B -spline  [14].  By  using  adaptive  tetrahedral  meshing 
and  a  mesh  untangling  technique,  an  algorithm  was  developed 
to  construct  a  trivariate  T-spline  representation  of  genus-zero 
solids  [4]. 

It  is  still  a  challenging  problem  to  automatically  create  poly¬ 
cubes  for  high  genus  geometry  and  use  them  in  constructing 
analysis-suitable  trivariate  T-splines.  In  this  paper,  we  utilize  a 
harmonic  field  defined  over  the  input  mesh  to  build  the  poly  cube 
automatically,  and  then  construct  the  rational  solid  T-spline  over 
the  polycube.  We  include  pillowing  and  quality  improvement 
techniques  to  guarantee  that  the  obtained  solid  T-spline  can  be 
used  for  analysis  directly. 

3.  Algorithm  Overview 

As  shown  in  Figure  1,  there  are  three  main  stages  for  con¬ 
structing  a  solid  T-spline  from  a  given  boundary  triangle  mesh 
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with  arbitrary  genus  topology.  From  the  input  triangle  mesh, 
we  first  compute  a  harmonic  scalar  field  defined  over  the  mesh, 
extract  the  geometry  topology,  and  then  generate  the  polycube 
with  the  same  topology.  Adopting  the  polycube  as  the  paramet¬ 
ric  domain,  we  build  the  valid  T-mesh  in  the  second  stage  and 
construct  solid  T- splines  in  the  last  stage. 


parametric  directions,  which  may  lie  on  an  edge  or  a  face.  In 
solid  T- splines,  an  edge  T-junction  is  a  T-junction  which  lies 
on  one  edge,  such  as  node  M  in  Figure  2(d).  A  face  T-junction 
is  a  T-junction  lying  on  one  face,  such  as  node  P  in  Figure  2(e). 


Fig.  2.  Four  types  of  nodes  in  solid  T-meshes.  (a)  Regular  node;  (b)  partial 
extraordinary  node;  (c)  extraordinary  node;  (d)  edge  T-junction;  and  (e)  face 
T-junction. 


Fig.  1.  An  overview  of  the  solid  T-spline  construction  algorithm  from  the 
given  boundary  triangle  mesh  with  arbitrary  topology. 

The  polycube  generation  stage  consists  of  two  steps: 

(i)  Harmonic  Field  Calculation  -  We  build  a  smooth  har¬ 
monic  scalar  field  defined  over  the  input  mesh.  Based  on 
this  field,  we  compute  its  gradient  field  and  an  orthogo¬ 
nal  vector  field; 

(ii)  Handing  Critical  Points  -  From  the  harmonic  field,  we 
determine  all  the  critical  points  where  the  first  partial 
derivatives  vanish.  These  critical  points  include  extreme 
points  and  saddle  points.  We  design  different  methods  to 
deal  with  each  type  of  point  in  order  to  build  the  poly  cube. 
The  polycube  edges  are  traced  along  the  gradient  and 
isocontour  directions. 

Based  on  the  poly  cube,  the  T-mesh  is  constructed.  There 
are  four  different  kinds  of  nodes  in  the  solid  T-mesh:  regular 
nodes,  partial  extraordinary  nodes,  extraordinary  nodes  and  T- 
junctions.  A  regular  node  is  a  node  around  which  the  local 
T-mesh  is  a  structured  mesh,  like  node  A  in  Figure  2(a).  Both 
partial  extraordinary  nodes  and  extraordinary  nodes  are  irreg¬ 
ular,  and  they  can  be  distinguished  using  reflection  edges.  Re¬ 
flection  edges  are  a  pair  of  adjacent  edges  with  one  common 
node,  and  the  set  formed  by  all  the  elements  sharing  one  edge 
is  topologically  symmetric  with  the  set  of  elements  sharing  the 
other.  For  example,  AB  and  AC  in  Figure  2(b)  are  a  pair  of 
reflection  edges.  A  partial  extraordinary  node  is  an  irregular 
node  about  which  some  but  not  all  of  its  adjacent  edges  have 
reflection  edges,  like  node  A  in  Figure  2(b).  An  extraordinary 
node  is  an  irregular  node  about  which  none  of  its  adjacent 
edges  has  a  reflection  edge,  such  as  node  A  in  Figure  2(c).  A 
T-junction  terminates  a  row  of  control  points  in  one  or  more 


The  T-mesh  construction  stage  consists  of  four  steps: 

(i)  Parametric  Mapping  -  We  build  a  parametric  mapping 
between  the  input  triangle  mesh  and  the  boundary  of  the 
generated  polycube; 

(ii)  Octree  Subdivision  and  Projection  -  An  initial  T-mesh  is 
obtained  by  an  octree  subdivision  of  the  polycube  and 
each  node  on  the  poly  cube  boundary  is  projected  onto 
the  boundary  surface  based  on  the  mapping; 

(iii)  Pillowing  and  Quality  Improvement  -  We  insert  one  pil¬ 
lowed  layer  on  the  boundary  and  improve  T-mesh  quality 
by  smoothing  and  optimization; 

(iv)  Handling  Extraordinary  Nodes  and  Partial  Extraordinary 
Nodes  -  In  order  to  obtain  a  gap-free  T-mesh,  we  apply 
templates  to  each  extraordinary  node  and  partial  extraor¬ 
dinary  node  in  the  initial  T-mesh. 

In  this  stage,  we  use  the  polycube  as  the  parametric  domain 
for  the  solid  T-spline  construction.  In  the  octree  subdivision 
step,  we  choose  the  existing  T-junction  parametric  values  to 
subdivide  each  octree  cell  as  much  as  possible,  instead  of  al¬ 
ways  using  the  central  parametric  value.  Based  on  the  valid  T- 
mesh,  the  knot  vectors  for  each  node  are  determined  by  travers¬ 
ing  T-mesh  faces  and  edges  [17],  and  the  solid  T-spline  is  con¬ 
structed  based  on  the  rational  T-spline  definition  [22].  Bezier 
elements  are  extracted  to  serve  as  the  primary  computational 
objects  in  isogeometric  analysis.  See  [2,  15]  for  elaboration. 

4.  Polycube  Generation 

In  this  section,  we  discuss  the  detailed  algorithm  of  the  poly¬ 
cube  generation  from  the  input  boundary  triangulation.  The 
polycube  must  be  constructed  in  a  geometrically  approximate 
and  topologically  equivalent  way.  To  achieve  this  goal,  we  first 
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compute  a  harmonic  function  defined  over  the  input  mesh  T, 
derive  two  orthogonal  fields  from  the  harmonic  function,  and 
extract  the  topology  structure  of  T  by  examining  critical  points. 

4.1.  Harmonic  Field  Calculation 

To  extract  the  topological  structure  of  the  input  mesh  T  c 
R3 ,  we  construct  a  harmonic  function  /  :  T  —>  R,  such  that 

A/  =  0  (1) 

subject  to  the  boundary  condition  that  vertices  in  the  predefined 
set  S min  and  Smax  have  the  minimum  and  maximum  values.  A 
is  the  Laplace  operator.  Smin  and  Smax  are  either  given  by  the 
user,  or  they  can  also  be  determined  by  selecting  the  bottom¬ 
most  and  top-most  points  on  T.  Basically,  computing  such  a 
harmonic  function  is  to  assign  a  scalar  value  to  each  vertex  in 
T.  For  a  triangle  mesh,  we  use  the  discretization  of  the  Laplace 
operator 

A  fi  =  L  Wj(f(.Vj)  -  /(Vi))  =  0,  (2) 

Wi 

where  Vu  Vj  e  T,  Wij  is  the  weight  and  A)  is  the  number  of 
vertices  adjacent  to  Vi.  Here,  we  choose  the  weights,  wtj  = 
cot  aij  +  cot  fiij,  where  or*/  and  fiij  are  the  opposite  angles  of  the 
edge  Vi  -  Vj.  From  the  discretization  of  the  Laplace  operator, 
the  scalar  function  /  can  be  obtained  by  solving  a  linear  system. 
Figure  3(a)  shows  the  scalar  function  defined  over  the  “Eight” 
model.  The  red  region  has  the  maximal  scalar  value  and  the 
blue  region  has  the  minimal  scalar  value. 

From  the  scalar  field  /,  we  compute  both  the  gradient  gi  =  V/ 
and  one  orthogonal  vector  field  g2  (the  isocontour  field).  Due  to 
the  properties  of  a  harmonic  scalar  field,  the  obtained  direction 
fields  are  guaranteed  to  be  smooth  and  free  of  extraneous  critical 
points.  For  one  triangle  (Vi,  Vj,  14),  suppose"^  is  the  unit  normal 
vector.  The  gradient  vector  g\  =  V/  is  obtained  by  solving  the 
following  linear  system  [3]: 


Vj-Vi 

fj  fi 

Vk-Vj 

bd  = 

fk-fj 

It 

0 

an  edge  AB  and  we  simply  follow  this  edge.  In  Figures  3(b-c), 
the  red  edges  are  the  newly  obtained  flow  line  segments.  The 
next  step  is  to  consider  vertex  B  and  trace  the  flow  line  from  it 
using  the  same  procedure. 


Fig.  3.  (a)  Harmonic  scalar  field  with  saddle  points  rendered  in  pink  for  the 
“Eight”  model,  (b-c)  Two  cases  for  tracing  the  gradient  flow  line  from  a 
given  vertex  A.  (b)  Regular  case;  and  (c)  edge  case. 


In  addition,  based  on  the  scalar  field,  we  can  then  determine 
all  the  critical  points  of  /,  that  is,  those  points  whose  partial 
derivatives  vanish.  These  points  include: 

-  Minimal  point  -  Points  in  set  Smin', 

-  Maximal  point  -  Points  in  set  Smax; 

-  Splitting  saddle  point  -  Points  where  the  geometry  splits; 

-  Merging  saddle  point  -  Points  where  the  geometry  merges. 


The  field  g2  is  along  the  isocontour  directions  of  /.  Hence,  for 
one  scalar  value,  we  simply  find  out  its  isocontour  to  obtain 
the  vector  field  g2  for  the  triangle  mesh.  Once  we  obtain  the 
two  orthogonal  vector  fields  g\  and  g2 ,  we  can  trace  along 
the  flow  lines.  A  flow  line  is  a  piecewise-linear  curve  defined 
over  the  mesh  whose  edges  are  along  one  of  the  vector  fields. 
There  are  two  cases  for  tracing  the  gradient  flow,  the  regular 
case  and  the  edge  case.  As  shown  in  Figures  3 (b-c),  vertex  A 
is  the  starting  point,  the  green  arrows  represent  the  gradient 
direction  for  each  incident  triangle,  and  we  wish  to  trace  the 
flow  line  from  A  by  walking  across  the  incident  triangles  along 
the  gradient  direction.  For  the  regular  case  in  (b),  starting  from 
A  we  extend  the  gradient  line  by  crossing  one  of  its  incident 
triangles.  In  this  case,  we  add  one  new  vertex  B  to  advance  the 
flow  line.  For  the  edge  case  in  (c),  the  flow  field  converges  on 


Fig.  4.  Configurations  around  a  regular  (a),  minimal  (b),  maximal  (c),  and 
saddle  point  (d). 

Figure  4  shows  the  configurations  of  a  regular,  minimal,  max¬ 
imal  and  saddle  point.  The  red  points  denote  the  vertices  with 
a  larger  scalar  value  compared  with  the  center  point  and  the 
blue  points  denote  the  vertices  with  a  smaller  scalar  value.  As 
shown  in  Figure  4(a),  around  one  regular  point,  vertices  on  one 
side  all  have  a  larger  scalar  value  and  vertices  on  the  other  side 
all  have  a  smaller  scalar  value.  For  a  minimal/maximal  point, 
all  the  vertices  surrounding  it  have  a  larger  or  smaller  scalar 
value.  For  a  saddle  point,  the  sign  changes  alternately  along  its 
circumferential  direction. 

Discussion:  The  harmonic  field  is  controlled  by  the  user- 
defined  minima  and  maxima  constraints,  it  then  affects  the  poly- 
cube  generation  and  alignment  of  the  solid  T- spline.  To  ensure 
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(a) 


(b) 


(c) 


(d) 


Fig.  5.  (a-b)  Handling  the  minimal  level;  and  (c-d)  handling  one  splitting  saddle  point  S . 


the  harmonic  field  follows  the  geometric  shape  of  the  input  sur¬ 
face,  generally  it  is  better  to  place  these  constraints  on  the  tips 
of  the  geometry  and  also  consider  symmetry.  For  example,  in 
the  “Eight”  model  (Figure  3(a)),  we  assign  the  constraints  at 
the  top  and  the  bottom. 

4.2.  Handing  Critical  Points 

We  need  to  handle  various  critical  points:  minimal,  maximal, 
splitting  saddle  and  merging  saddle  points.  We  first  compute 
all  the  scalar  levels  or  isocontours  at  which  there  are  critical 
points.  Let  denote  the  isovalue  of  Li  and  (Dj  denote  the 
corresponding  isocontour.  For  level  Li  with  saddle  points,  two 
sets  of  isocontours,  CL  and  Ct,  are  computed  by  using  the 
isovalue  with  a  small  perturbation  5.  For  the  minimal  or 
maximal  level,  only  one  isocontour  is  computed. 

Suppose  level  Li  contains  one  minimal  point.  At  level  Li , 
four  seed  points  (P\pj  =  0, ...,3)  are  chosen  for  each  closed 
curve  as  shown  in  Figures  5 (a-b),  which  correspond  to  the  four 
lower  comers  of  the  cube  Ci  with  the  unit  parametric  length. 
The  parametric  value  of  the  other  vertices  lying  on  this  isocon¬ 
tour  will  be  computed  using  the  chord-length  parameterization. 
From  these  seed  points  P\  .  we  trace  the  gradient  flow  line  un¬ 
til  it  intersects  with  the  isocontour  (CL  {  of  the  next  level  L/+ 1 
at  P“j  (Li+ 1  may  contain  a  saddle  or  maximal  point).  The  red 
curves  in  Figure  5  are  traced  gradient  lines,  and  the  black  ones 
are  isocontours.  The  four  vertices  PL-  -  (J  =  0, . . . ,  3)  correspond 
to  the  four  upper  corners  of  Ci.  The  traced  four  curves  are  then 
mapped  onto  the  four  vertical  edges  of  the  cube,  and  P\  .  and 
Put  j  serve  as  the  eight  corners.  Then  the  poly  cube  construction 
process  advances  to  the  next  level. 

For  level  Li  with  a  splitting  saddle  point  S  as  shown  in 
Figures  5(c-d),  we 

(i)  Parameterize  the  lower  isocontour  (CL  using  Pl.  .  =  . 

as  four  corners  (assuming  the  associated  cube  before 
splitting  is  C;_ i); 

(ii)  Find  the  shortest  path  from  the  splitting  saddle  point  S 
to  (CL,  get  the  intersection  node  Qo  on  (CL,  and  calculate 
point  Qi  on  the  opposite  edge  with  the  same  parametric 
value; 

(iii)  Compute  the  shortest  path  between  Qo  and  Q\  (see  the 
blue  curve  in  Figure  5(c),  the  path  can  not  contain  any 
edge  on  this  isocontour); 


(iv)  Determine  two  sets  of  points  on  (CL ,  Qq  -  Q J  and  <2q  -  Q\ , 
which  have  the  shortest  distance  from  Qo~Q\  to  (CL ; 

(v)  Construct  two  cubes  and  C]  by  using  Pli0-Qo-Q\-Pli3 
and  Qo~P\  r^j2-Gi  as  the  lower  comers,  respectively; 

(vi)  continue  tracing  the  gradient  flow  until  the  flow  line  in¬ 
tersects  the  isocontour  at  the  next  level. 

Basically,  for  one  splitting  saddle  point,  we  aim  to  find  one 
isoparametric  line  to  split  one  cube  into  two.  Here,  the  isopara¬ 
metric  line  connecting  Qo~Q\  is  used  to  split  the  upper  face 
of  the  cube  Ci- 1  as  shown  in  Figures  5 (c-d).  Similarly,  for  each 
merging  saddle  point,  we  use  the  same  procedure  to  find  one 
isoparametric  line  to  merge  two  neighboring  cubes  and  ensure 
they  match  with  each  other  seamlessly. 


(a)  (b)  (c) 

Fig.  6.  Polycube  construction  and  mapping  results  for  the  “Eight”  model,  (a) 
The  constructed  polycube  in  the  parametric  space;  (b)  the  poly  cube  in  the 
physical  space;  and  (c)  the  polycube  mapping  result. 

By  dealing  with  all  the  critical  points,  the  poly  cube  is  con¬ 
structed  level  by  level.  We  always  map  the  isocontour  onto 
the  horizontal  isoparametric  edges  of  the  polycube  and  map 
the  gradient  flow  lines  onto  the  vertical  edges.  Figure  6  shows 
one  poly  cube  construction  result  for  the  “Eight”  model.  The 
red  lines  in  Figure  6(b)  denote  the  curves  corresponding  to  the 
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edges  of  the  poly  cube.  This  polycube  construction  algorithm 
does  not  consider  the  symmetry  property  of  the  input  geom¬ 
etry.  However,  if  Smin  and  Smax  are  given  symmetrically,  the 
obtained  polycube  and  the  constructed  solid  T-spline  can  be 
symmetric  for  a  symmetric  input  geometry. 

Discussion:  Here,  we  only  consider  Morse  saddle  points 
whose  multiplicity  is  1.  Morse  saddle  points  are  handled  by 
splitting  one  cube  into  two,  or  merge  two  cubes  into  one.  In 
general,  for  saddle  points  of  any  multiplicity,  one  cube  may  be 
split  into  an  arbitrary  number  of  cubes.  If  there  are  two  or  more 
saddle  points  on  one  level,  multiple  sets  of  <2o-<2i  need  to  be 
computed. 

5.  T-mesh  Construction 

The  T-mesh  construction  stage  aims  to  build  one  valid  T- 
mesh  from  the  given  boundary  triangulation  and  the  constructed 
polycube.  There  are  four  main  steps  in  this  stage:  parametric 
mapping,  adaptive  octree  subdivision  and  projection,  pillowing 
and  quality  improvement,  handling  extraordinary  and  partial 
extraordinary  nodes. 

5.1.  Parametric  Mapping 

This  step  aims  to  build  a  one-to-one  parametric  mapping 
between  the  input  triangle  mesh  T  and  the  boundary  of  the 
obtained  polycube,  which  serves  as  the  parametric  domain  of 
the  solid  T-spline.  From  the  constructed  polycube,  we  have  the 
correspondence  between  the  traced  gradient  or  isocontour  lines 
and  the  edges  of  the  polycube.  Based  on  the  traced  lines  on 
the  triangle  mesh,  we  can  divide  the  input  mesh  into  N  sub¬ 
meshes,  Tl  (i  =  0, . .  . ,  A),  where  N  is  the  number  of  boundary 
faces  of  the  polycube.  Each  sub-mesh  is  associated  with  one 
face  of  the  poly  cube,  Tc\i  =  0, . .  . ,  A).  We  then  use  the  surface 
parameterization  to  map  each  sub-mesh  Tl  to  its  corresponding 
poly  cube  face  Tc  • 

For  each  sub-mesh,  we  first  map  its  boundary  to  the  boundary 
of  Tc  by  a  chord  length  parameterization.  The  parameteriza¬ 
tion  for  the  interior  vertices  is  calculated  by  solving  a  linear  sys¬ 
tem  formed  by  the  harmonic  equation  ^  w;7(/( Vj)  -  /(VO)  = 

Wi 

0,  where  Wij  =  cot  or^-  +  cot fiij.  For  the  curve  shared  by  two  ad¬ 
jacent  sub-meshes,  we  use  the  same  parameterization.  Figure 
6(c)  shows  the  mapping  result  for  the  “Eight”  model.  To  guar¬ 
antee  a  conformal  boundary  between  two  neighboring  cubes, 
we  always  choose  the  same  parameterization  for  all  the  edges 
shared  by  them.  Note  that  the  two  neighboring  cubes  A  and  B 
do  not  share  faces.  There  are  two  duplicated  faces  in  the  para¬ 
metric  space  but  they  are  separate  in  the  physical  space. 

5.2.  Adaptive  Octree  Subdivision  and  Projection 

An  initial  T-mesh  is  generated  by  applying  an  adaptive  octree 
subdivision  to  the  poly  cube  C  and  projecting  to  the  boundary. 
For  each  cube  in  C,  we  create  one  hexahedral  element,  using 
the  same  parametric  length  and  considering  the  physical  length 
difference  in  three  directions.  Then  we  obtain  a  root  T-mesh  for 


the  whole  polycube  and  treat  it  as  one  single  piece,  instead  of 
treating  each  cube  separately.  Starting  from  the  root  T-mesh,  we 
subdivide  one  element  into  eight  smaller  ones  recursively  to  get 
the  refined  initial  T-mesh  after  projection.  For  each  boundary 
element,  we  check  the  local  distance  from  the  T-mesh  boundary 
to  the  input  triangular  mesh,  and  subdivide  the  element  if  the 
distance  is  greater  than  a  given  threshold  s. 

Each  obtained  node  has  both  parametric  and  physical  coor¬ 
dinates.  The  parametric  coordinates  represent  its  position  in  C. 
For  each  boundary  node,  the  physical  coordinates  are  its  asso¬ 
ciated  position  on  the  triangular  mesh  based  on  the  polycube 
mapping.  The  physical  coordinates  of  each  interior  node  are  cal¬ 
culated  by  a  linearly  interpolation.  Note  that  the  three  isopara¬ 
metric  planes  are  not  always  in  the  middle.  If  one  element  con¬ 
tains  T-junctions,  the  parametric  values  of  the  T-junctions  are 
used  to  subdivide  the  element.  If  there  are  more  than  two  T- 
junction  parametric  values  in  one  direction,  the  one  closest  to 
the  middle  is  used.  For  example,  the  purple  element  in  Figure  7 
has  one  T-junction  on  the  £-edge  and  two  T-junctions  on  the  77- 
edge.  For  the  £  direction,  the  parametric  value  £1  is  used,  while 
for  the  77  direction  772  is  used  during  subdivision,  because  it  is 
closer  to  the  parametric  middle.  The  dash  lines  are  inserted  to 
refine  this  element.  In  this  way  we  can  minimize  the  number 
of  T-junctions  in  the  initial  T-mesh. 


% 

c 

% 

- 1 - 

1 

1 

1 

1 

1_  -J 

) 

_ A _ 

L. , 5  * 


Fig.  7.  T-mesh  subdivision.  Green  nodes  denote  T-junctions. 

The  octree  subdivision  and  projection  processes  continue  un¬ 
til  the  local  distance  from  each  boundary  element  to  the  input 
mesh  is  less  than  s  and  each  element  has  at  most  one  edge  T- 
junction  for  an  edge,  or  one  face  T-junction  for  a  face.  Figure 
8  gives  one  result,  (a)  shows  the  parametric  coordinates  for  the 
nodes  in  the  T-mesh  and  (b)  shows  their  physical  coordinates. 

5.3.  Pillowing  and  Quality  Improvement 

To  improve  the  T-mesh  quality,  we  adopt  the  pillowing, 
smoothing  and  optimization  techniques.  Pillowing  is  a  sheet- 
insertion  technique  that  inserts  one  layer  around  the  boundary 
[23].  Here  we  insert  one  pillowed  layer  for  the  initial  T-mesh, 
which  helps  to  improve  the  T-mesh  quality  and  the  T-spline 
surface  continuity. 

Figure  9  illustrates  the  pillowing  operation  for  the  poly  cube. 
The  cubes  C;  are  rendered  in  different  colors.  As  mentioned 
earlier,  C\  and  C2  do  not  share  a  face  in  the  parametric  space. 
In  (b),  we  use  two  separate  faces.  The  dark  blue  lines  show  the 
pillowed  layer.  In  pillowing,  each  boundary  face  is  duplicated 
to  form  one  pillowed  element,  and  each  pillowed  element  has  at 
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(a)  (b) 


(d) 


(a)  (b) 


Fig.  9.  The  polycube  before  (a)  and  after  (b)  pillowing, 
where  Nj  is  a  trilinear  shape  function.  The  scaled  Jacobian  is 


Fig.  8.  The  subdivision,  projection,  pillowing  and  optimization  results  for  the 
“Eight”  model,  (a)  The  subdivision  result  in  the  parametric  domain;  (b)  the 
initial  T-mesh;  (c-d)  Interior  of  the  T-mesh  before  and  after  pillowing  and 
optimization. 


most  one  face  lying  on  the  boundary.  For  the  pillowed  layer,  the 
edge  knot  interval  is  a  predefined  constant  along  the  pillowing 
direction,  which  stays  the  same  for  the  other  two  directions. 

As  shown  in  Figure  9(b),  after  pillowing  the  yellow  cor¬ 
ner  nodes  become  interior  extraordinary  nodes.  The  nodes  on 
the  blue  polycube  edges  become  interior  partial  extraordinary 
nodes.  The  red  corner  nodes,  pillowed  from  the  yellow  nodes, 
are  extraordinary  nodes  on  the  boundary  surface.  The  con¬ 
structed  T-spline  surface  is  C° -continuous  around  these  red 
boundary  extraordinary  nodes  up  to  the  2-ring  neighborhood, 
C1 -continuous  from  the  2-ring  to  3-ring  neighborhood,  and  C2- 
continuous  everywhere  else.  Note  that  after  pillowing,  the  sur¬ 
face  continuity  across  the  polycube  edges  is  improved  from  C° 
to  C2.  The  interior  region  is  C°-continuous  around  each  interior 
extraordinary  node  until  the  3 -ring  neighborhood.  For  the  inte¬ 
rior  region  across  the  poly  cube  edges,  the  continuity  is  C°  until 
the  2-ring  neighborhood,  C1  from  2-ring  to  3 -ring  neighbor¬ 
hood,  and  C2  everywhere  else.  Since  here  we  introduce  some 
extraordinary  and  partial  extraordinary  nodes,  we  use  a  local 
parameterization  for  each  element  in  the  following  steps. 

After  pillowing,  smoothing  and  optimization  are  used  to  im¬ 
prove  the  T-mesh  quality.  For  smoothing,  each  node  is  moved 
toward  its  mass  center,  and  for  optimization  each  node  is  moved 
toward  an  optimal  position  that  maximizes  the  worst  Jacobian. 
For  one  T-mesh  element,  the  Jacobian  is  defined  as 


J  =  det(Jiw)  — 
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where  Jm(‘,  1)  and  Jm( -,2)  represent  the  first,  second 

and  last  column  of  the  Jacobian  matrix,  Jm ,  respectively.  To 
handle  T-junctions  during  smoothing  and  optimization,  we  add 
some  “virtual  nodes”  to  refine  the  local  region  and  convert  the 
local  T-mesh  to  a  hexahedral  mesh.  Figure  8(d)  shows  the  result 
after  pillowing  and  optimization  for  the  “Eight”  model. 

5.4.  Handling  extraordinary  and  partial  extraordinary  nodes 

Unlike  regular  nodes,  extraordinary  nodes  or  partial  extraor¬ 
dinary  nodes  may  introduce  gaps  to  the  solid  T-spline.  In  this 
step,  we  apply  templates  given  in  [22]  to  make  the  T-mesh  gap- 
free.  Figure  10(a)  shows  the  template  for  a  partial  extraordinary 
node,  in  which  the  magenta  edge  has  a  reflection  edge.  Figure 
10(b)  shows  one  general  template  for  an  extraordinary  node. 


Fig.  10.  The  general  template  for  a  partial  extraordinary  node  (a)  and  an 
extraordinary  node  (b).  The  magenta  node  is  a  partial  extraordinary  node, 
the  red  one  is  an  extraordinary  node,  and  the  blue  ones  are  inserted  nodes. 
The  magenta  edge  is  the  edge  with  a  reflection  edge  and  the  red  edges  have 
zero  knot  interval. 


6.  Solid  T-spline  Construction  and  Bezier  Extraction 
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In  this  stage,  we  aim  to  infer  the  local  knot  vectors  for  each 
node,  build  the  rational  solid  T-spline  from  the  T-mesh,  and 
then  extract  embedded  Bezier  elements  [2,  15].  The  concept 


of  rational  T-splines  was  given  in  [22],  with  basis  functions 
satisfying  a  partition  of  unity  by  definition.  The  rational  solid 
T-spline  is  defined  as 


;=o 


(g,ri,OeCl, 


(6) 


;=o 


where 


RM,r],0  = 


Nf(^N^T)Njiq 

rj=o^)N’’(r])N(j(0 


(7) 


is  the  rational  B-spline  basis  function,  vf,  V  and  Aif  are  B- 
spline  basis  functions  defined  by  the  local  knot  vectors  at  node 
Ci  which,  for  degree  d  =  3,  are  given  by  £  =  [£0,  £i,  fe,  fo, 
^74  L  ^li  =  [^7/0  5  ^7/ 1 5  ^2,  ^3,  ^7/4  ]  and  ^  =  [^0?  (Tz‘l>  £z'2?  £z'3?  (Tz'4  ]  • 

From  the  T-mesh,  we  first  compute  the  knot  intervals  for 
each  node  by  traversing  T-mesh  edges  and  faces.  Note  that  we 
repeat  knots  whenever  we  meet  an  extraordinary  node  during 
the  traverse.  Then  for  each  domain,  we  determine  all  the  nodes 
with  non-zero  basis  functions  and  use  them  to  build  the  solid 
T-spline  element.  The  entire  solid  T-spline  model  is  built  by 
looping  over  all  the  local  domains. 

Bezier  extraction  provides  a  finite  element  representation  of 
T-splines  for  isogeometric  analysis.  Similar  to  [23],  we  first 
compute  the  “Bezier  mesh”,  which  includes  all  the  reduced 
continuity  lines  by  adding  all  the  knot  vector  inference  lines 
or  isoparametric  planes  for  T-junctions  and  L-j unctions  to  the 
T-mesh.  Then  for  each  element  in  the  Bezier  mesh,  the  trans¬ 
formation  matrix  Me  between  the  T-spline  basis  functions  and 
the  Bezier  basis  functions  is  calculated  satisfying: 


Bet  =  MeBeh , 


(8) 


where  Bet  is  the  vector  formed  by  the  nonzero  T-spline  basis 
functions  in  this  element  and  Beb  is  the  vector  formed  by  the 
trivariate  Bezier  basis  functions. 


7.  Results  and  Isogeometric  Analysis 

We  have  applied  the  construction  algorithm  to  several  models 
(Figures  11-14).  The  constructed  solid  T-spline  is  tricubic  and 
C2 -continuous  except  in  the  vicinity  of  partial  extraordinary 
and  extraordinary  nodes.  Statistics  for  all  the  tested  models  are 
shown  in  Table  1.  The  Bezier  Jacobian  is  calculated  using  the 
scaled  Jacobian  at  the  eight  Gauss  quadrature  points  for  each 
Bezier  element.  The  algorithm  is  efficient  and  all  the  results 
were  computed  on  a  PC  equipped  with  an  Intel  X3470  processor 
and  8  GB  main  memory. 

The  Isis  model  has  genus  zero  and  its  polycube  only  con¬ 
tains  one  single  cube.  For  the  kitten  model,  there  are  three  sad¬ 
dle  points  (one  merging  saddle  point  and  two  splitting  saddle 
points).  We  only  consider  two  of  them,  and  the  splitting  sad¬ 
dle  point  close  to  the  maximum  point  is  skipped  to  get  a  more 
simplified  polycube.  The  sculpture  model  has  genus  two  and 
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(a)  (b)  (c) 

Fig.  11.  The  solid  T-spline  construction  result  for  the  “Eight”  model,  (a)  The 
constructed  solid  T-spline  and  T-mesh;  (b)  the  extracted  solid  Bezier  mesh 
with  some  elements  removed  to  show  the  interior  (blue)  and  one  pillowed 
layer  (magenta);  and  (c)  the  isogeometric  analysis  result. 

five  saddle  points.  Again,  we  skipped  the  one  near  the  max¬ 
imal  level.  The  time  used  for  each  model  not  only  depends 
on  the  input  mesh  size,  but  also  depends  on  the  T-mesh  size, 
which  is  determined  by  the  given  surface  error  tolerance  and 
the  complexity  of  the  topology  and  geometry.  One  advantage 
of  this  harmonic  function  based  T-spline  construction  is  that 
the  isoparametric  lines  are  basically  aligned  with  the  geometric 
structure  of  the  model. 

We  have  developed  a  3D  isogeometric  analysis  solver  for 
static  mechanics  analysis,  which  uses  rational  T-splines  as  the 
basis,  and  we  used  it  to  test  the  obtained  solid  T-spline  models. 
For  all  the  models,  we  fix  all  the  control  points  on  the  bottom 
and  apply  unit  displacement  on  the  top.  The  Young’s  modulus 
E  =  HAGPa,  and  the  Poisson’s  ratio  v  =  0.3.  The  obtained 
displacement  results  are  given  in  Figures  1 1(c),  12(e),  13(g)  and 
14(g).  From  these  results,  we  can  conclude  that  the  obtained 
rational  T-splines  can  be  used  in  isogeometric  analysis  directly 
and  reasonable  simulation  results  can  be  obtained. 

8.  Conclusions 

We  have  presented  a  new  algorithm  to  construct  solid  T- 
splines  for  arbitrary-genus  geometries  from  boundary  triangu¬ 
lations.  Our  method  is  efficient  and  the  resulting  solid  T-spline 
is  analysis-suitable  with  C2 -continuity  except  for  the  local  re¬ 
gion  around  a  few  irregular  nodes.  In  this  paper,  we  only  con¬ 
sider  the  geometries  with  Morse  saddle  points,  and  as  part  of 
our  future  work  we  will  extend  the  algorithm  to  geometries 
with  arbitrary  saddle  points.  We  will  also  consider  engineering 
designs  with  sharp  features. 
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Fig.  12.  Isis  model  with  genus  zero,  (a)  The  input  boundary  triangle  mesh;  (b)  the  harmonic  scalar  field;  (c)  the  constructed  solid  T-spline  and  T-mesh;  (d) 
the  extracted  solid  Bezier  elements;  (e)  the  isogeometric  analysis  result;  (f)  the  constructed  polycube  and  the  mapping  result;  (g)  the  subdivision  result  for  the 
parametric  domain;  and  (h)  the  extracted  solid  Bezier  mesh  with  some  elements  removed  to  show  the  interior  (blue)  and  one  pillowed  layer  (magenta). 
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Table  1.  Statistics  of  all  the  tested  models 


Model 

Genus 

Input  mesh 

(vertices,  elements) 

T-mesh  Extraordinary  nodes 

nodes  (Boundary,  Interior) 

Interior  partial 

extraordinary  nodes 

Bezier 

elements 

Bezier  Jacobian 

(worst,  best) 

Time 

(s) 

Isis 

0 

(5,863,  11,722) 

9,310 

(8,  8) 

244 

5,335 

(0.12,  1.0) 

52.0 

Kitten 

1 

(5,377,  10,754) 

4,825 

(12,  12) 

140 

2,883 

(0.08,  1.0) 

14.8 

Eight 

2 

(766,  1,536) 

5,735 

(16,  16) 

200 

1,440 

(0.10,  1.00) 

8.5 

Sculpture 

2 

(8,635,  17,276) 

10,549 

(16,  16) 

252 

7,072 

(0.09,  1.00) 

41.5 

Fig.  13.  Kitten  model  with  genus  one.  (a)  The  input  boundary  triangle  mesh;  (b)  the  harmonic  scalar  field;  (c)  the  constructed  polycube  and  the  mapping  result; 
(d)  the  subdivision  result  for  the  parametric  domain;  (e)  the  constructed  solid  T-spline  and  T-mesh;  (f)  the  extracted  solid  Bezier  elements;  (g)  the  extracted 
solid  Bezier  mesh  with  some  elements  removed  to  show  the  interior  (blue)  and  one  pillowed  layer  (magenta);  and  (h)  the  isogeometric  analysis  result. 
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Fig.  14.  Sculpture  model  with  genus  two.  (a)  The  input  boundary  triangle  mesh;  (b)  the  harmonic  scalar  field;  (c)  the  constructed  polycube  and  the  mapping 
result;  (d)  the  subdivision  result  for  the  parametric  domain;  (e)  the  constructed  solid  T-spline  and  T-mesh;  (f)  the  extracted  solid  Bezier  elements;  (g)  the 
extracted  solid  Bezier  mesh  with  some  elements  removed  to  show  the  interior  (blue)  and  one  pillowed  layer  (magenta);  and  (h)  the  isogeometric  analysis  result. 
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SINTER  The  “Eight”  model  is  from  the  AIM  @  SHAPE  Shape 

Repository. 
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